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Abstract 

In this paper we present recent developments concerning a Cell-Centered 
Arbitrary Lagrangian Eulerian (CCALE) strategy using the Moment Of 
Fluid (MOF) interface reconstruction for the numerical simulation of multi- 
material compressible fluid flows on general unstructured grids in cylindri- 
cal geometries. Especially, our attention is focused here on the following 
points. First, we propose a new formulation of the scheme used during the 
Lagrangian phase in the particular case of axisymmetric geometries. Then, 
the MOF method is considered for multi-interface reconstruction in cylin- 
drical geometry. Subsequently, a method devoted to the rezoning of polar 
meshes is detailed. Finally, a generalization of the hybrid remapping to cylin- 
drical geometries is presented. These explorations arc validated by mean of 
several test cases that clearly illustrate the robustness and accuracy of the 
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Figure 1: Multi-material CCALE algorithm flowchart. 



In this work, we consider the simulation of multi-material compressible 
flows on unstructured meshes in cylindrical geometry. For this, we adopt 
an ALE description [T3] that has the great advantage to combine the best 
features of both Eulerian and Lagrangian approaches. Indeed, this choice is 
not only well adapted to naturally track free surfaces and interfaces between 
different fluids as purely Lagrangian methods, but also to handle flow distor- 
tion as Eulerian methods. Here, a CCALE [TUl [TT] approach is particularly 
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considered whose the main elements are as follow. 

As depicted on figure Fig|T} the first step of the algorithm relies on an ex- 
plicit Lagrangian phase in which the physical variables and grid are updated 
thanks to a slightly modified version of the Explicite Unstructured Cell- 
Centered Lagrangian HYDrodynamics (EUCCLHYD) scheme [I9l EHl El] in 
cylindrical coordinates. Recently, new investigations have been made about 
cell-centered Lagrangian schemes [SI [7] . The scheme presented in this paper 
is a modified version of the area weighted finite volume scheme of jlH]. Then, 
multi-material fiows treatment is done thanks to specific interface capturing 
method. This choice allows to track the volume fraction of each material used 
for the thermodynamical closure relying on the equal strain rates assump- 
tion. This approach is quite simple to implement and to use and remains 
sufficient in almost cases [HI [23]. This, leads to constant evolution of the 
volume fraction during the Lagrangian phase. Such an approach allows to 
reconstruct with accuracy the interface between each material. In this con- 
text, many development have been done for 2D Cartesian geometries. First, a 
previous version of the CCALE algorithm solving two-material compressible 
fiows using a Volume Of Fluid (VOF) have been proposed in [HI E]. Then 
an extension to Moment Of Fluid (MOF) approach has been considered to 
enhance multi- material (more than two components) fiows in [HI |T0] . Subse- 
quently, a rezoning phase is realized. It consists in moving the Lagrangian 
nodes to improve the geometric quality of the grid [U]. Finally, the physi- 
cal variables are conservatively interpolated from the Lagrangian grid onto 
the new rezoned one during the remapping phase. Here an extension of the 
hybrib remapping [1] to cylindrical geometries is introduced. We want to 
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notice that in ALE framework using cell-centered formulation, this phase is 
straightforward. In the lines of these works, the main goal of this paper is 
to extend the CCALE-MOF algorithm to treat both Cartesian and cylindri- 
cal geometry. To this end, several modifications are given to the algorithm 
previously presented. In a first part, we propose a new formulation of the 
numerical scheme introduced in [19] for treating axisymmetric geometries 
during the Lagrangian phase. To build this scheme, an area-weighted formu- 
lation of the Lagrangian system of equations is proposed. Then, this system 
of equations is discretized using a cell-centered finite volume (FV) scheme. 
Contrary to [19] in which fluxes are directly deduced from the Geometric 
Conservation Law (GCL) constraint, here a simpler formulation that gives 
similar results is retained. These two main choices lead to a robust first- 
order scheme conservative for the total energy that has the great advantage 
to preserve spherical symmetry for one- dimensional flow on uniform angular 
polar grids. The high order extension has been performed using the Gener- 
alized Riemann Problem (GRP) described in [19]. To treat interface flows, a 
MOF interface reconstruction method is retained in the sequel. Once again, 
the difficulty here is to propose a natural and consistent adaptation of this 
approach able to treat axisymmetric interface flows. To this end, formula- 
tions of the moments needed to track interface are revisited for cylindrical 
coordinates as in |2]. This leads to an accurate and second order interface 
reconstruction method that allows to treat multi-material (more than two) 
interfaces in the lines of [8]. The third part of this study is dedicated to 
recent enhancement of the rezoning algorithm to improve the mesh quality 
during computation especially on polar meshes. As it is done in [TIH [TT| . 
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mesh rezoning is based on the Condition Number Smoothing (CNS) [TB] al- 
gorithm on unstructured meshes. Moreover, when used for polar meshes, it 
is well known that CNS algorithm pushes the nodes toward the origin dete- 
riorating the mesh quality. To avoid this drawback, the main idea developed 
in this paper is to adapt CNS algorithm to polar grids. Then, extension to 
unstructured grids (Cartesian-polar) is also explored. Finally, a generaliza- 
tion of the remapping procedure to cylindrical geometries is proposed. Here, 
an efficient method adapted to multi-material flows is presented. The main 
idea is to use an hybrid remapping that combine the main advantages of 
the swept-face and multi-material cell- intersection remapping as in |H ITO] . 
Finally, a specific attention is done to polynomial integration that preserves 
the method efficiency. 

The paper is structured as follows. We detail in the second section a 
new formulation of the first-order area weighted Lagrangian scheme used for 
axisymmetric geometries. Further extensions to high-order are notably de- 
tailed in p!9]. Afterwards, the extension of the MOF axisymmetric interface 
reconstruction method is presented for treating multimaterial fiows. Then, 
we describe the General Condition Number Smoothing (GCNS) algorithm 
for unstructured meshes. Finally, the description of the new hybrid remap- 
ping procedure for cylindrical geometry is done. For a complete description 
of the CCALE-MOF method see [TUl [H], except new advances presented 
in this paper. Then presentation of numerical experiments is made in Sec- 
tion 4. They demonstrate not only the robustness and the accuracy of the 
present methodology but also its ability to handle successfully complex two- 
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dimensional multi-material fluid flows notably computed for axisymmetric 
geometries. Finally concluding remarks and perspectives about future works 
are given in the last section. 

2. Lagrangian phase in axisymmetric geometry 

In this part, an extension of the cell-centered Lagrangian scheme pUf fIT\ 
is presented for the numerical simulation of compressible flows in pseudo- 
Cartesian geometries for unstructured meshes as in pjj. This choice has 
the great advantage to treat both axisymmetric and Cartesian geometries. 
In this paper, a new and simple formulation of the scheme introduced in 
[19] for flrst-order approximation is proposed. To this end, an area weighted 
formulation of classical Lagrangian equations is flrst introduced. Then these 
equations are discretized with a node-centered approximate Riemann solver. 

2.1. Governing equations 

During the Lagrangian phase, the rates of change of volume, mass, mo- 
mentum and total energy are computed assuming that discretized volumes 
move following the flow. Thus, each arbitrary volume V{t) depending on the 
time t > moves satisfying the following system of equations 




(1) 




(2) 
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where — is the Lagrangian derivative and p, U, P, E are respectively the 
at 

density, velocity, pressure and total energy. In addition, this system is closed 
thanks to an equation of state (EOS) as 

P = p{p,£), (5) 

with the internal energy e defined as e = E — |Up/2. At last, we have local 
kinematic equation 

U, X(0)=X°, (6) 



~dt 



with X the location of a point of the control volume surface S{t), at time 
t > and X*^ its initial value. This equation is equivalent to (|2| also known 
as geometric conservation of law (GCL). 



2.2. Area-weighted formulation 
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Figure 2: Notations related to the pseudo-cartesian grid. 

For defining the differential operators used in the system of Lagrangian 
equation ([l])-(|4]) a pseudo-Cartesian reference frame {0,X, y} for the or- 
thonormal basis (ex,ey) is used (see FigJ2]). Thus each point is localized 



by means of its positions X and TZiY) = 1 — a + aY the pseudo-radius. 
When a = 0, the Lagrangian equations for Cartesian geometry are recov- 
ered, otherwise for a = 1 this corresponds to axisymmetric equations. In 
this way, axisymmetric geometry is obtained from Cartesian one through a 
rotational symmetry about the X-axis. This imphes that the volume V{t) is 
generated by the rotation of the area A{t) about the X-axis. In consequence, 
the element volume dV writes as dV = TZdA with dA = dXdY the element 
area in the pseudo-Cartesian frame. In the same manner, the control surface 
S{t) delimiting V{t) is obtained through the rotation of L{t) the boundary 
of A{t) and the surface element is given by dS = TZdL. Note that we have 
omitted the 27i factor in the evaluation of the element volume. 

In a such framework, the velocity divergence and the pressure gradient 
read as follows 



Using the previous definitions and after some calculations using the Green's 
formula, it is possible to rewrite ([l])-(|4]) at least in two different ways. The 
first one, obtained without any approximation is the control volume formu- 
lation. When discretized this formulation leads to a conservative scheme for 
both equations of energy and momentum, and satisfies the local semi-discrete 
entropy inequality. However, as shown in [TP] it does not preserve symme- 
tries. Consequently, an area-weighted formulation is adopted here leading 




where U* = (m, v) 



(7) 



and 
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to a conservative scheme for energy equation that respect spherical geome- 
tries. This formulation is deduced from the control volume one assuming 
that momentum equation ^ is written in Cartesian geometry. Like this, the 
area-weighted formulation for the Lagrangian equations reads 

m— {V)+n / PNrfL = 0, (10) 

m—{E)+ / PN -niJdL = 0, (11) 

where m = Jy pdV represents the mass of the volume V . Each physical 
variable per unit of mass {E, U) is noted as 0, and has its mass density 
mean value defined by (0) = ^ /y p(pdV. The average 71 corresponds to 
ratio 71 = ^. In such case, as m = pV, the momentum equation is solved in 
Cartesian geometry. For Cartesian case V = A, we recover TZ = 1. Further 
details on the derivation of this system are available in [19]. 

2.3. Numerical scheme 

Thereafter, we recall briefly the first order cell-centered Lagrangian scheme 
introduced in [TH]. To this goal, similar notations as [1011121 [22] are employed 
in the sequel. Let us consider a set {i^cjceN of non-overlapping polygonal cells 
that approximates A{t). Each cell noted Qc is assigned a single index c. Each 
vertex of the cell c is labeled with the index p and is localized thanks to its 
coordinates Xp = (Xp,!^)* in the pseudo-Cartesian frame. In addition, we 
introduce V{c) the list of the vertices belonging to the cell Qc and C{p) the list 
of the cells sharing the vertex p. These two sets are counterclockwise ordered. 
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Let us introduce p~ and p'^ the previous and the next nodes with respect to 
p in V{c). We denote by L~^,Lp^ the half length of the edges [pp~], [pp~^]- 
Similar notations are used for the normals outward N^^ and N^^. Finally, 
the corner normal Lp^Npc is given by Lp^Np^ = L+^N^^ + L~^N~^. All these 
notations have been displayed in FlGjSj 




Lpc^pc — LpJ^p^ + LpJ^p^ 



= Vc/Ac 



i\,p 1 p 



Figure 3: Notations for the cell-centered scheme. 



The first order spatial approximation of (|9|-(11) is obtained considering 
local integrals on each cell Qc rotated about the X-axis. The mass rric of the 
cell Qc is rric = pdV and each flow variable (as total energy, velocity) 
is averaged over each cell through the formula 



m, 



p(j)dV, 
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named cell-centered value. Then, we have 

m,jV, + n, J2 Fpc = 0, (12) 

Yl Fp,•7^pUp = 0. (13) 

In addition, the mesh is moved through the local kinematic equation given 
at each node by 

^ = Up for t > and Xp(0) = Xj, (14) 

with Up and respectively the velocity and the position of a node p at 
initial time. In the previous equations, Fp^ is the numerical flux at each 
node p of each cell c deflned by 

Fpe = Lp,P,Npc-Mp,(Up-Uj, (15) 

with Up the velocity at the point p and Pc the mean value of the pressure in 
the cell c. The 2x2 matrices Mpc and Mp are deflned as 

Mp, = Z, (L- N- ® N- + L+ N+ ® N+ ) , and Mp = J] Mp,. (16) 

cec(p) 

Where, we introduce the "swept mass flux" [9] associated to the isentropic 
sound speed a, that is 

Zc = Pcttc- (17) 

This is nothing but the acoustic impedance. As it has been demonstrated in 
[22] the total energy and momentum conservation is equivalent to 

Fpc = 0. (18) 

cec(p) 
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Finally using (15), the nodal velocity Up is deduced from (18) by solving 
the linear system 

MpUp = ^ (VP.Npe + Mp,Ue). (19) 

cec(p) 

In [19], the numerical fluxes used in the discretization of (|9| and (11) are 



chosen for satisfying the local GCL constraint (14). Here, we rather adopt a 



more simple approach that give similar results. Since ( 14 ) is explicitly solved 
for moving the mesh in time, there is no need to solve Thus, each cell 
volume Vc is directly deduced from (14). Thereby, it is possible to choose for 



the numerical flux in (11) a simple form as in (13) with TZp = 1 — a + aYp. 

Concerning, the momentum equation the mean value TZc is equal to the 

— K 
discrete ratio TZc = —r- 

Ac 

Let us note that this new formulation of the area-weighted discretiza- 
tion relies on a node-centered solver which is exactly the same as the one 
developed in [20] for two-dimensional Cartesian geometry. However, the 
present spatial discretization does not satisfy rigorously the GCL compat- 
ibility requirement. In what follows, we will assess the discrepancy of our 
discretization to the GCL by analyzing the corresponding discrete divergence 
operator. The discrete divergence operator that corresponds to the present 
scheme writes as 



{^■U)c=^Y. ^Pi^pc^pc + L^eN^c) ■ Up, 



(20) 



where TZp denotes the pseudo-radius of vertex p. It is shown in [22], |T9] that 
the discrete divergence operator deduced from the GCL reads 

iy.UtcL^^ J2 l[(27^p + 7^p)L^-Ar + (27^p + 7^+)L+A^+].C/p. (21) 



Vc ^ 3 
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If the time evolution of the position vector, Xp, of vertex p is governed by 
the trajectory equation [14], then one can prove that the time rate of change 



of the cell volume, V^, satisfies 



^ "^^^ - (V ■ C/)^^^. 



K dt 



Subtracting (20) and (21) leads to 



(v-c/)e-(v-C7),^^^ = ^ J2 [(7^p-7^,)L^-A^-,+(7^,-7^+)L+A^+].c/p. 

(22) 

Knowing that the summation in the previous equation is cyclic, shifting the 
index in the second term of the right hand-side yields 

(V ■ t/), - (V ■ t/)^^^ = ^ E l(K - T^p)LtcN;^] ■ (Up. - Up). (23) 

p(iV{c) 

In case of a one-dimensional spherical flow on an equi-angular polar grid, 
the right-hand side of the previous equation is equal to zero. To prove this 
result, let us consider a quadrangular cell of an equi-angular polar grid. The 
proof proceeds in the following two steps: 

• Either p and p^ are located on the same angular sector and thus the 
nodal velocity Up and Up+ are colinear to the direction of the angular 
sector which is orthogonal to the unit outward normal AT^. Hence, 
(C7,+ - Up) ■ N+ = 0. 

• Or p and p"*" are located on the same cercle of radius R, then the 
Cartesian components of their nodal velocities reads as 

/cos A (cosi9 + M) 

Up = U{R) , C7p+ = U{R) 

ysinej ysm{9 + A9) 
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Here, 9 denotes the angle of the angular sector, U{r) is the module of 
the one-dimensional velocity field, and /\9 is size of the angular sector. 
A straigthforward computation shows that 



we obtain that {Up+ — Up) ■ AT^ = 0. 

This ends the proof. This result shows that our new area-weighted dis- 
cretization satisfies rigoroulsy the GCL compatibility requirement for one- 
dimensional spherical flows on equi-angular polar grids. 

3. MOF multi-material interface reconstruction phase in axisym- 
metric geometry 

The method used in this work to reconstruct interfaces, is the MOF ap- 
proach well adapted for treating multi- materials interface problems [U |8]. 
Indeed, such a method enables to capture more accurately interfaces than 
the classical VOF strategy and allows the treatment of general multi-material 
flows (more than two materials) [ini [E] • This method has been recently ex- 
tended in cylindrical geometries, for a single interface problem [2]. Here, 
extension to multi-material interface reconstruction phase to cylindrical co- 
ordinates is considered. 




Knowing that the unit outward normal is given by 
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3.1. Moment of fluid method 

The main idea of MOF is to track each fluid in a cell using the zeroth 
and first moments [8]. Given these two moments, interface is linearly re- 
constructed insuring volume conservation. To this end, interface update is 
done minimizing the discrepancy between the given moments and the recon- 
structed moments of the polygon behind the interface. One should note that 
no information from neighboring cells is required. This method is exact for 
linear interfaces and is second order accurate for smoothly curved ones. In 
the context of multi-material configurations, one has to face to material or- 
dering when reconstructing interface. The method presented here, allows to 
automatically determine the order of materials by constructing all the pos- 
sible combination and choosing the sequence that leads to the configuration 
where the reconstructed moments are the closest to the given ones. The main 
difference between cylindrical and planar geometry relies in the definition of 
the different moments. Since the interface reconstruction is done under vol- 
ume conservative assumption, the zeroth moment M^^ of the /c-th fluid in 
each cell c is obviously given by 



with the cell volume K = /j^ TZdA. 

Contrary to the zeroth moment, the first moment can be defined without 
any specific requirement. Thus, it is possible to compute them in the two 




(24) 



from this moment we can deduce the volume fraction 




(25) 
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following different manners. In the one hand we can use the natural extension 
to axisymmetric geometries 

Ml^ = I n^dA, (26) 
and from this moment we deduce the pseudo-centroid 

ML , ^ 

Xfc,c = 77^, with \4,c = Kttfc.c- (27) 

This pseudo-centroid for a matter of simplicity will be called here the ax- 
isymmetric centroid. 

On the other hand it can also be done with a planar definition as follows 

<f = / XrfA, (28) 



k.c 



and thus planar centroid will be obtain from 

Xt=^, (29) 

^k,c 

where A^ ^ is the area of the k-th fluid in the cell c. 

Since this interface reconstruction method is coupled to our Lagrangian 
hydrodynamics scheme it requires to update the volume fractions and ma- 
terial centroids. Using the equal strain assumption, the volume fractions do 
not evolve during the Lagrangian step (see [TT] for more details). However, 
the centroid locations are given from the Lagrangian step using a barycentric 
combination of the new positions of the mesh nodes as done in [TU] . 

3.2. Numerical validation 

The main goal of this section is to compare the results given by both 
axisymmetric and planar formulations of the centroids on several static test 
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cases in one cell. As in [5], we consider three different mixed-cell layouts 
that are filament (without junction), T-junction and Y-junction. The first 
two configurations correspond to C^-serial partitions whereas the third is not. 
In the considered test cases, the parameter x corresponds to the radius of 
the circles defining the interfaces. Two values are considered with x = 1 ^'^^ 
X = 64. In addition, the computation domain is reduce to the cell [0; 1] x [0; 1] 
(see figures FlG]4]and Fig]5]). 

In the first case, with x = 1, we notice small differences for the filament 
case, no notable difference on the T- Junction but the Y-junction results 
for axisymmetric and planar formulations present distinct interface positions 
due to a different ordering of the materials. For a large radius x = 64, the 
curves are reduced to piecewise linear interfaces. Then, the result using both 
formulations are very close to each other. For the two first cases filament 
and T-junction, the results are exact. Regarding the Y-junction, it remains 
a good approximation. These results illustrate the capability of both planar 
and axisymmetric centroid formulation for MOF to treat accurately multi- 
material problem. Nevertheless, for consistency with the global cylindrical 
coordinate formulation, the axisymmetric formulation for the centroids is 
retained in the sequel. 
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Filament T-junction Y-j unction 




Figure 4: MOF interface reconstruction test for three materials. From the top to the 
bottom: the true partitions for x — 1 ^-nd their MOF reconstructions obtained with 
planar and axisymmetric centroids. From the left to the right: filament, T-junction and 
Y-junction configurations. 
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Filament T-junction Y-junction 




Figure 5: MOF interface reconstruction test for three materials. From the top to the 
bottom: the true partitions for x = 64 and their MOF reconstructions obtained with 
planar and axisymmetric centroids. From the left to the right: filament, T-junction and 
Y-junction configurations. 
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4. Rezoning phase improvement for polar meshes 

The rezoning phase introduced in [TOl [1^ consists in moving the La- 
grangian grid to improve its geometric quahty. The objective of this part is 
to extend this approach to polar meshes. To this end, the proposed proce- 
dure relies on two main steps. The first phase is dedicated to compute the 
smoothed grid from the Lagrangian one through CNS method. Then the 
final mesh is deduced from the smoothed one by a relaxation procedure to 
keep the rezoned grid as close as possible to the Lagrangian grid in order 
to insure computation accuracy and avoid unphysical mesh rezoning. In the 
sequel one should note that rezoning is formulated only for planar geometry 
in the frame {0, X, Y}. 

For the sake of readability, in the rest of the paper the quantities without 
any accent a are associated to Lagrangian mesh. After the rezoning step 
we use a*"^^, and finally after relaxation the quantities related to the rezoned 
mesh are noted with the tilde accent a. 

4-1. General condition number smoothing (GCNS) 

As it is pointed out in the introduction, CNS approach is well adapted 
to rezone Cartesian meshes but it still suffers from drawbacks for polar ones. 
Indeed, in this case the mesh seems to collapse (like an implosion) to the ori- 
gin. To circumvent this difficulty, it has been proposed to modify the CNS 
algorithm using specific weight associated to the mesh geometry [25] that 
controls mesh rezoning with regards to the radius for example. Nevertheless, 
this approach is not completely satisfactory. First, it strongly depends on 
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the choice of the weight, that may affect the quahty of the mesh which can 
be shifted in the opposite direction to the origin for example. Furthermore, 
there is still a residual compression near the origin due to singularity at this 
point. In conclusion, it does not preserve a uniform polar mesh. For this 
reason, a different strategy is presented here. The main idea developed here 
is to apply the CNS rezoning algorithm in (r, ^)-coordinate system. In fact, a 
polar mesh initially expressed using a Cartesian coordinates {X, Y) leads to a 
structured Cartesian mesh in (r, ^)-coordinates. Here, a general presentation 
of the algorithm is made for unstructured meshes. 



Assuming that the resulting mesh from the Lagrangian phase is unfolded 
(otherwise untangling procedure is used to correct invalid cells |26]) . Thus, 
the proposed algorithm consists for polar meshes in three different steps as 



depicted on Fig. 17 For the sake of simplicity, we consider in the sequel only 
the case of Cartesian and polar structured meshes. 




T 



CNS 



rezoning 



T 



-1 




lag nlag\ 



(r 



rez nrez\ 



Figure 6: Rule representation for GCNS algorithm for a polar mesh. 



The first step, is dedicated to the mapping between Cartesian and polar 
coordinates. To this end, consider c a given cell of the Lagrangian grid for 
(X, F)-coordinates, p E V{c) a node of this cell. Notation used in the sequel 
are depicted on Fig|7} The mapping between a point p G c of Cartesian 
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p' 




T = R, 



P 



Xp = (X,y) Xp = (r,0) 

• mapped nodes, o extrapolated nodes 
Figure 7: Notations and mapping between cartesian and polar coordinates. 

coordinates Xp = (X, YY to Xp = (rp, 6pY in polar ones is done through the 
following linear transformation 



Xp — TXp, 



where TfX^ 



def 



h - l^h + /3R(Xp) with 13 e {0, 1}. For /3 



(30) 
1 then 



TfX. 



R(Xp) with the rotation matrix 



RfX. 



cos 



sm 



sm 



cos (6'p) 



using the definition 9p = arctan '^p ~ V-^p ~^ ^v' ^^^^ 

/3 = 0, this formula leads to Cartesian rezoning with the identical transfor- 
mation Xp = Xp. When mapping (X, F) to (r, 6*), the origin node has to 



be specifically treated. Indeed the transformation (30) is not defined for this 



point. Then as it is needed in the rezoning algorithm in the (r, ff) frame, the 
origin node is defined by a mapping of the first row on r = axis (see Fig|7]). 
Note that these nodes are not used for the final backward mapping. 

The second step is the GCNS algorithm. It is based on a minimization 
problem of a local functional that controls the quality of the mesh. As done 
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in [inilll], one has to distinguish boundary nodes and internal node for which 
the smoothing procedure is different. 

For internal nodes, let us introduce as in ^2] the condition number for (r, 6)- 
coordinates that writes 

1 1 V 1 12 I 1 1 V 1 12 

«:(Jcp) = , (31) 

where Xpp± = Xp — Xp±, and A^p = det(Jcp) is the area of the triangle 
delimited by {p,p^,p^} in the rezoned grid and Jcp = [Xpp+| — Xpp-] the 
2x2 Jacobian matrix associated to each corner at vertex p of cell c. Thanks 
to this condition number we define the local function associated to the node p 

Fp(Xp) = «^(Xp), (32) 
cec(p) 

Finally, the new position Xp^^ is obtained by the minimization of the local 
function Fp using the first step of a Newton algorithm. This leads to the 
formula 

X;- = Xp - H-i(Xp)VFp(Xp), (33) 

where H;^^ and VFp are respectively the Cartesian 2x2 Hessian matrix and 
gradient related to the local functional Fp. 

For boundary nodes, the rezoned position Xp*^^ of p is computed in consistent 
way with the GCNS algorithm. To this end, Xp*^^ is given thanks to a second- 
order interpolation Bezier curve [11] leading to 



X;^^ = Xp{s''') = (1 - {s'''f)Xp- + 2(1 - s'-^^)Xi + (s^^^)'Xp+, (34) 

where Xj such that A:'p(l/2) = Xp. Furthermore, the parameter s^^^ is com- 
puted to minimize Fp{Xp{s)) (for more details on this procedure see [H]). 
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Finally, the third step consists in backward mapping between X^*^^ and 



Xp*^^ using (30), where the inverse of the transformation matrix is taken 
equal to T-i(Xp) =^ h -Ph + P [R(Xp)]~^ with /3 G {0, 1} and [R(Xp)]"^ 
the inverse rotation matrix . 

4-2. Relaxation algorithm 

The relaxation algorithm consists in making a convex combination be- 
tween rezoned grid obtained from GCNS step and its location after La- 
grangian step. This reads for each mesh node p by: 

Xp = Xp + UpiX'/' - Xp), with Up G [0, 1], 

where Xp is the new mesh node position after the complete rezoning phase. 
The coefficient Up is computed as a function of the right Cauchy-Green tensor 
associated to the Lagrange grid deformation over a time step (for details see 

mm)- 

4-3. Numerical validation 

In this section, we compare results obtained by the GCNS algorithm 
to those obtained for classical CNS for the rezoning of uniform polar and 
unstructured meshes. 

Uniform mesh. First, we consider an uniform polar mesh made of 20 x 10 
elements see FlG|8]-(a). Results obtained after 100 iterations for the clas- 
sical and general smoothing are presented on FlGjSj For each method the 
relaxation coefficient Up is taken equal to 1. As already mentioned, the clas- 
sical smoothing does not converge on polar mesh and implies the collapse 
of cell layers to the origins (see FlG|8}(b)). However, for the GCNS, the 
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result obtained (see FlG.|8]-(c)) is converged. The mesh initially uniform, is 
not modified at the end of the computation. This clearly illustrates the good 
behavior of our smoothing algorithm. 




(a) (b) (c) 

Figure 8: Smoothing of a static polar grid 16 x 10: (a) initial grid; Smoothed grids after 
100 iterations: (b) CNS, (c) GCNS. 

Unstructured mesh. Now, rezoning for an unstructured mesh is studied. 
Let us consider a mesh made of 175 quadrangular cells as depicted on Fig|9} 
(a). When applying the full Cartesian rezoning to the mesh, similar obser- 
vations as previously can be made. It suffers from an implosion of central 
cells to the origin and does not converge (see FlG]9]-(b)). For the full GCNS 
algorithm, one can see after convergence, the formation of mesh distortion 
on the square region and a polar mesh far from the center (see Fig|9]-(c)). 
Nevertheless, it is possible to improve this rezoning. Thus, the main idea 
developed in the sequel is to apply the GCNS rezoning algorithm differently 
for a node belonging initially to a Cartesian or polar region of the mesh. 
To this end, the transformation T between {X, Y) and (r, 6) coordinates is 
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(a) (b) (c) 

Figure 9: Smoothing of a static unstructured grid: (a) initial grid; Smoothed grids after 
100 iterations (b) CNS, (c) GCNS after 100 iterations. 

modified in the following way 

T(Xp) I, - /5(Xp)/, + /3(Xp)R(Xp), (35) 

with 

1 if xo G 

/3(Xp) 



p 

if x° e 



p 

where P™*^ and P^"' are the sets of nodes that belong to the Cartesian and 
respectively polar region of the mesh at the initial time. These regions are 



represented thanks to red and blue color (see FlG.lO-(a)) for the considered 



mesh. Nodes localized at the frontier between the polar and Cartesian meshes 



(black nodes on Fig. 10 -(a)) can be considered either polar, or Cartesian. 
As represented on FlG|l0|-(b,c), both possibilities are tested. The obtained 
results illustrate that the Cartesian choice remains better contrary to the 
polar one that introduce mesh distortion. 
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(a) (b) (c) 

Figure 10: Smoothing of a static unstructured grid: (a) initial grid with Cartesian (blue) 
and polar (red) rezoning regions; Smoothed grids after 100 iterations (b) GCNS with 
interfacial polar rezoning, (c) GCNS with interfacial Cartesian rezoning. 

5. Hybrid remapping in axisymmetric geometry 

During the remapping phase, the physical unknowns (density, velocity, 
total energy) computed thanks to the Lagrangian step are conservatively 
remapped from the Lagrangian mesh to the rezoned one. To this end, an ex- 
tension of the Hybrid Remapping Algorithm for multi-material flows [101115111] 
to cylindrical geometry is proposed here. This strategy consists in the follow- 
ing two steps. First a swept-faced remapping is used to treat cells and nodes 
localized far from the interface. Then, a cell-intersection-based method [TT] 
is applied to the cells and nodes in the neighborhood of the interface. In 
this way, this approach combines the ability of the cell-intersection method 
to remap the interface and the efficiency of the swept flux approach for the 
other cells that significantly reduce the global computing cost of the method. 
As done previously, in the perspective of general use of the method, a global 
formulation including both Cartesian and axisymmetric framework is pre- 
sented. 
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We assume in the sequel, that there is no topology change of the mesh, 
the cells of the Lagrangian and rezoned grids are respectively designed by Qc 
and Qc- 

5.1. Multi-material cell-intersection-based (MCIB) remapping 
Y 

a \ 

I- - - - 



ey 




ex 



X 



Figure 11: Notations for MCIB method. 

The main goal of remapping is as follows. Given the piecewise constant 
representation of the physical variables per unit of volume (p, pU, pE) noted 
i^c = Pc<Pc in each cell of the Lagrangian grid, we want to compute its equiv- 
alent ipc in each cell of the rezoned grid given as 

1 



%jjc = — I p(f)TZdA, 



(36) 



with Vc the volume of the cell Vt^. Contrary to single fluid approach, here 
the rezoned values ipc can not be computed directly in each cell c. In fact, 
one has to take into account multi-material aspects. 



28 



First of all, let us introduce some notations. Each material of the flow 
noted k occupies the polygon Qk,c C Qc, within the MOF framework, such 
that Qc = [_J^k,c and is characterized by its partial mass, density, pressure, 

k 

internal energy and variables per unit of volume (total energy, momentum) 
whose averaged values in each sub-cell are respectively rrik^c, Pk,c, Pk,c, £k,c and 
'4'k,c = Pk,c4>k,c with (pk,c the partial velocity or energy per unit of mass. 

Thus, for multi-material flow, the main idea of remapping is not to di- 
rectly compute the global rezoned quantities ipc but the partial rezoned ones 
noted ■ilJk,c- This is particularly true for the MCIB method that is dedicated 
to treat cell in the interface neighborhood. To this end, we first propose 
a second order reconstruction \E'fc,c(X) of ipk^c over each Lagrangian cell c 
through the piecewise linear function 

^fc,c(X) = V^fc,, + (V^fc),(X - Xfc,e), (37) 

where (V^fc)c denotes the constant gradient of '^k,c within cell c computed 
thanks to a least-squares approach. Finally X^^c is the centroid related to 
the k-th fluid in the cell c given by 

Xfc,, = [ nXdA. (38) 

Thanks to these notations, the remapped value for MCIB is given by 

^k,c=^ J2 f ~ '^^k,cdA, (39) 

where the intersection polygons Qk,d H Qc are computed thanks to a specific 
triangulation of the mesh. The procedure is detailed in [TT]. The set C(c) con- 
tains the cells including c that share at least one node with the cell c. At last, 
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the partial volume defined on the rezoned cell .s V>,. = E / ^ ^'^A. 

In the context of MOF reconstruction, one has to define additional quanti- 
ties as the partial remapped mass corresponding to material k. It is computed 
as fhk^c = 'f>k,cOik,cVk,c with the volume fraction 



(40) 



thus the partial volume can be also expressed as V^^c = K«A;,c- In addition, 
each material centroid position is defined thanks to 



5.2. Pure cell swept-face (PCSF) remapping 



(41) 




k,c 



k,c+ 



Figure 12: Notations for swept face-based method. 



As explained before, the PCSF remapping is used only to treat single 
fiuid cells. In this context, one should remark that fi^ = ^k,ci thus the mean 
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value ipck is given through 

V'fc,c = V'fc,c+ I n^kjdA, (42) 

with Af the quadrangular signed area swept by the face / of a cell c between 
the Lagrangian grid and the rezoned grid delimited by the ordered nodes of 



coordinates {Xp, Xp, Xp+, Xp+} (refer to fig Fig. 12). We note J^(c) the set 



of the faces / of a cell c. In addition, '^kj is the upwind value given by 

{^kc+ if v4f > 
' (43) 
'^k,c otherwise. 

with the neighbor cell of c through the face /. During this step the volume 
fractions ak,c = o:k,c do not change as we consider single fluid cells and the 
material centroid can be updated directly from the geometry X^^c = Xc 
where X^ is the centroid of the cell Qc- 

5.3. Integration strategy 

For both PCSF and MCIB remapping, one has to compute several surface 
integrals, on polygons where the integrand is a polynomial function of (X, Y). 
This can be done using a triangulation of these areas. Nevertheless, this is 
expensive. Here, we rather adopt a more efficient method as in [23]. In this 
context, integrals are simplified using Taylor decomposition of the polynomial 
integrand and Green's formula leading to compute circular integrals over the 
edges of the polygons defining the integration areas. For further details on 
integral computations see 
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Figure 13: Hybrid remapping principle in one-dimension case. 

5.4- Hybrid remapping algorithm 

In this part, we detail the hybrid remapping algorithm that is summarized 



on Fig. 13 To this end, let us introduce and N^'^ the sets of nodes and 
in the same manner and C^^ the sets of cells respectively used for PCSF 
and MCIB remapping. Here A/"^^ collects mixed nodes belonging to cells 



that contain the interface or are on this interface (white nodes on Fig. 13) 



despite contains the pure ones (black nodes on Fig. 13). In addition, 
is the set of mixed cell that include cells intersected by the interfaces and 
their neighbors by nodes. Finally, contains the cells that have at least 
one node in . 

The hybrid remapping procedure consists in performing the following steps. 
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1. PCSF step. In this step we first move the pure nodes included in 
and we remap the quantities in cell c belonging to . Thus, we 
have ^fc,c = (Pfc,c, P-Efc,c5 pUfc J using relation Q and m^^c, 5fc,c, ^k,c 



have il)k,c = {Pk,c, pEk,ci P^k,c) using relation Q and m^^c, ak,c, '^k,c 
for each cell c E . 
2. MCIB step. Now, the mixed nodes in Af^ are moved and the ilJk,c = 
{pk^c, pEk^c, pV f^c) remapped thanks to (39) and mk^c,(ik,c,^k,c are 
computed for cells c G C*^. 

Since C^^ fl 7^ {0}, one should note that cell included in this intersection 
remapped at each step of the algorithm. 



are 



At the end of remapping, only the partial values of the physical variables 
per unit of volume are known. A this step, a first point is to compute the 
physical variables per unit of mass. The remapped partial total energy is 
given using Ek^c = {pE)i:Jpk,c- However, this is different for the remapped 
partial velocity lJk,c- Indeed, as explained in the second part of this paper, 
the Lagrangian computation of the velocity is done in Cartesian geometry. 
For this reason, the remapped velocity is deduced from the {p\J)f^^ through 
Ufe,c = (pU)^ c/Pkc using the jo/anar remapped density and momentum given 



through (39) and (42) with 71 = 1. The second point is dedicated to the 
reconstruction of the global values required for the next Lagrangian step. To 
this end, a classical procedure is to use specific averages 

0c = i rhk,c4>k,c, (44) 
nic ^-^ 

k 

with the global mass and density deduced from 

fhk,cak,c and Pc = ^ Pk,cOik,c- (45) 



k k 
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At last, thermodynamical variables as pressure P and internal energy e are 
obtained thanks to specific thermodynamical closures as done in [II] . 

6. Numerical results 

We present in this section several numerical test cases performed using 
the CCALE-MOF computing procedure detailed in [101 E] including the 
various development proposed in this paper. In the sequel, all the materials 
are governed by perfect gas equation of state p = pe(7 — 1), where 7 stands 
for the polytropic index of gas. 

6.1. Sedov problem 




Figure 14: Initial grid and material positions for the Sedov problem. 

We present in this first section a Sedov problem for a point blast in a uni- 
form medium with spherical symmetry. We use this test case to compare our 
new formulation with the original EUCCLHYD scheme in pure Lagrangian 
and coupled to the CCALE-MOF procedure. The initial conditions are given 
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by {p^, P^,\J^) = (1,10~^,0) in a spherical domain of radius 1.2 except in 
the cell at the origin (0, 0) where an initial delta-function energy source is 
set through the pressure 

Par- = (7 - l)Por:^, 
Vor 

with Vor the volume of the origin cell and So = 0.851072 is the total amount 

7 

of released energy. The fluid has its polytropic index 7 equal to -. Contrary 

5 

to the original single material test case, we add here three artificial inter- 
faces, to test our multi-material CCALE-MOF algorithm. These interfaces 



are initially located for a radius equals to 0.1, 0.2 and 0.3 (see Fig. 14) 



Here we consider both Lagrangian and ALE computations for an initial un- 



structured mesh depicted on Fig. 14 This grid is obtained after one rezoning 
step, with Up = 1 of an unstructured mesh initially paved with 500 quad- 
rangular cells. Numerical results are depicted on Fig|T5] andFlGjT6] for a 
final time of tend = 1 and compared to the analytical solution computed us- 
ing self-similar arguments as done in [llj. It consists of a diverging shock 



wave whose front is exactly localized at radius R = 1. As it is illustrated 



on Fig. 15, the pure Lagrangian solutions are in good agreement with the 
analytical one for both approaches. We can notice that the new formula- 
tion is less dissipative as we reach a higher density level in the shock region. 
Indeed for the Lagrangian method as for the CCALE-MOF one the shock 



location is well resolved without any spurious oscillation (Fig. 16). In addi- 
tion, this simple problem underlines the robustness (better mesh quality near 
the origin) and accuracy (shock location) of the axisymmetric CCALE-MOF 
approach especially when considering multi-material flows whose interfaces 



are well captured thanks to the MOF reconstruction (see Fig. 16). 
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New Lagrangian scheme Original EUCCLHYD scheme 




Figure 15: Sedov problem. From the top to the bottom: Interface positions, density maps, 
density profiles defined as a function of the cell center radius compared to the analytical 
solution at final time step for the pure Lagrangian computation using new scheme (left) 
and the original EUCCLHYD scheme (right). 
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New CCALE-MOF EUCCLHYD CCALE-MOF 




Figure 16: Sedov problem. From the top to the bottom: Interface positions, density maps, 
density profiles defined as a function of the cell center radius compared to the analytical 
solution at final time step for the new CCALE-MOF procedure (left) and the EUCCLHYD 
CCALE-MOF procedure (right). 
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We point out that during the Lagrangian computation, non-convex cells 
appeared. This may lead to interface reconstruction failure when consider- 
ing multi-material flows. As illustrated by the previous numerical results, 
the proposed CCALE-MOF algorithm remains adapted to treat such config- 
uration without any difficulty demonstrating once again its robustness. 

6.2. Axi- symmetric triple point problem 



3 

(P2,P2,72,U2) '// 

=(0.125,0.1,1.5,0) % 



i Pi=l 



^ 71=1-5 
I Ui=0 



(P3,P3,73,U3) ^ 
, =(1.,0.1,1.4,0) I 

1 7 ex 



Figure 17: Axi-symmetric triple point problem : geometry and initial data. 

We consider in this part a three-material problem that corresponds to a 
three-state Riemann problem in an axisymmetric geometry. This problem 
has been wisely studied in Cartesian geometry and here we propose new re- 
sults for cylindrical geometry. The computational domain is rectangular and 
composed of three regions (blue, green, red) whose dimensions are depicted on 



Fig. 17, The top, left and right boundaries are closed thanks to walls. A sym- 
metry condition is applied to the bottom boundary corresponding to the X- 
axis axi-symmetry. Initially, the blue region contains a fluid with high pres- 
sure and density taken equal to (pi,Pi) = (1, 1). The green region contains 
a low density and pressure fluid whose initial state is (p2,P2) = (0.1,0.125). 
The third fluid in the red region, initially has a low pressure and an high 
density equal to (p3,P3) = (1,0.1). At the beginning of the computation, all 
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fluids are supposed to be at rest then Ui = U2 = U3 = 0. The blue and 
green material have the same polytropic index 71 = 72 = 1.5, despite the red 
one has 73 = 1.4. 




Figure 19: Axi-symmetric triple point problem. Mesh and material positions at t = 5 for 
ALE computation. 

The computation using the presented axisymmetric extension of the CCALE- 
MOF algorithm is made on a grid initially paved with 140 x 60 square cells 
until a final time tj = 5. For this simulation, comparison with a full La- 
grangian computation can not be performed since its suffers from important 
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mesh tangling as shown in [TS]. However comparison to full Eulerian simula- 
tions is done. In this case, nodes are moved to their initial positions during 
the rezoning step. Numerical results for both ALE and Eulerian methods 



representing interfaces and meshes are depicted on Fig. 18 19 As expected, 



since there is a shock wave with high speed that propagates from the heavy 
material (blue) to the light one (red), the interface is sheared at the triple 
point producing a Kelvin-Helmholtz like instability. Here, comparison to 
planar 2D computations (TD] demonstrates that axisymmetric geometry par- 
ticularly affects the vortex shape that is 3D. Although the global behavior of 
the solutions is very similar comparing ALE approach to the Eulerian one. 



6.3. Spherical Air-Helium shock/bubble interaction test 

Piston 



0.0445 



Ai 



/ 0.32 

(pi,Pi) = (0.182,105) 




0.65 



ex 



(P2,P2) = (1,10-5) 



Figure 20: Air-Helium shock/bubble interaction. Initial geometry and data. 

We deal in this part with the numerical simulation of the experiment of 
concerning the impact of a Mach 1.25 shock travelling through the air 
onto a spherical bubble of Helium. To this goal, let us consider a rectan- 
gular domain of dimensions [0, 0.65] x [0, 0.0445] initially full of Air of data 
(pi,Pi) = (0.182,10^) except in an half disc centered in (0,032) of radius 
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Fluid 


Polytropic index 7 


Molar mass ^ 


Air 


1.4 


28.963 


Helium 


1.648 


5.269 X 10-3 



Table 1: Air-Helium shock/bubble interaction: EOS parameters. 
0.0225 that contains Helium characterized by {p2,P2) = (1, 10^) as depicted 



on Fig. 20 Here, spherical geometry is obtain thanks to a rotation around the 
X-axis. EOS parameters for each fluids are stated on Tab|TJ Wall boundary 
and symmetry conditions are respectively chosen for the left, top boundaries. 
Despite, a piston-like condition is imposed to the right one for an incoming 
velocity equal to U* = ('U*,0). Here, the horizontal velocity u* is computed 
thanks to Rankine-Hugoniot conditions and is given by u* = —140.312 cor- 
responding to an incident shock moving at the velocity Dc = —467.707. 
The domain is initially paved with a structured cartesian grid composed of 
520 X 72 cells. Here, the bubble is directly initialized through the volume 
fraction on this mesh. Computations are done for the multi-material axisym- 
metric CCALE-MOF for a final time chosen equal totf = tj-|-600x 10^^ where 
ti = 657.463 X 10^^ corresponds to the time of the shock/bubble interaction. 
Here once again, simulations can not be achieved using pure Lagrangian 
framework due to the apparition of important mesh distortion. Numerical 
results associated to the Schlieren density profiles [13] and interface positions 



deduced from the MOF reconstruction are respectively depicted on Fig. 22 



and FlGj23| Let us note that each pictures are obtained thanks to an axial 
symmetry with respect to the X-axis. Comparisons between the Schlieren 
density profiles and the sadow-graphs of the experiment show a good agree- 
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ment, especially when observing the bubble shape deformations. Moreover, 
waves generated by the initial shock are well localized and illustrate multiple 
reflections and refractions especially on the bubble and the domain bound- 
aries. These main points clearly demonstrate the accuracy and the robustness 
of the method and validate the axisymmetric CCALE-MOF approach when 
computing spherical test-cases coming from experiment. 

6.4- Spherical implosion 



i 




Figure 21: Multi-mode implosion in spherical geometry. Initial geometry and data. 

The last test-case of this paper deals with the numerical computation of a 
spherical implosion as initially treated in [27] . The interest of this simulation 
is twofold. First, this is a realistic problem quite close to those encountered 
in Ignition Conflnement Fusion (ICF) simulation. Then, it allows to test 
the capability of the multi-material CCALE-MOF algorithm with hybrid 
rezoning. 
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t = ti + U5 X 10 





0.2 0.225 0.25 0.275 0.3 0.325 0.35 



t = ti + 223 X 10 




0.2 0.225 0.25 0.275 0.3 0.325 0.35 



t = ti + 350 X 10 





0.15 0.175 0.2 0.225 0.25 0.275 0.3 



t = ti + 600 X 10~^ 

Figure 22: Spherical Air-Helium shock/bubble interaction. Schlieren diagram of density. 
Axi-symmetric CCALE-MOF results (on the left) compared to experimental results (on 
the right) 12]. 
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t = ti + 20 X 10' 
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0.22 0.26 0.3 0.34 
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Figure 23: Spherical Air-Helium shock/bubble interaction. Mesh and material interface 
evolution after the shock hits the bubble at time ti — 657.463 x 10~^. 
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Here we focus on the treatment of perturbed interfaces where compressible 
Rayleigh- Taylor instabilities occur. 

Let us consider a spherical ball of light fluid (r G [0, 10]) initially surrounded 



by a shell of heavy fluid {R G [10, 12]) as depicted on Fig, 21 For both fluid 

5 

the polytropic index is the same 7; = 7/i = -• The initial pressures and 
densities are {pi,pi) = (0.05,0.1) and {ph,Ph) = (1,0.1). The implosion is 
driven by imposing the following pressure law on the dense shell boundary 



fit) 



10 iftG [0,0.5], 
12 -4t iftG [0.5,3]. 

Finally, the interface between the light and the heavy fluids is initially per- 
turbed according to the law 

r7 = rp(l + aoI^(rp)PKcos(0p)) 

with the damping factor 



r. - ri 



if rp G [ri,re 



1-^^^^ if, e[0,r,]. 

where r^^^ denotes the perturbed radius and oq is the amplitude of the per- 
turbation. Finally, Pi is the l^^ Legendre polynomial. In the sequel / = 10 
and several values of are considered from the non-perturbed case Oq = 0, 
to weakly and strongly perturbed one with respectively = 2 x 10~^ and 
ao = 1 X 10-^ 

Computations are made for two different meshes until the final time tf = 3. 
The first one is a polar grid displayed on FlGj24|-(left) composed of 90 x 40 
cells. Size of cells in the radial direction have been chosen respecting a mass 
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radial spacing deduced from the equivalent one-dimensional test case. The 
other grid, is obtained after an hybrid regularization for Up = 1 of an un- 
structured mesh initially paved with 3200 quadrangular cells respecting the 



mass radial spacing (see Fig. 24 -(right)). 



Non-perturbed case with = 0. As a first study, we test the behavior of our 
algorithm in axisymmetric geometries in pure Lagrange computation for both 



meshes. As shown on Fig, 25 numerical results for both meshes are similar. 



Nevertheless, the method remains faster on the unstructured mesh. Indeed, 
it has the advantage to not impose a drastic time step for computation due 
to triangular cells with high aspect ratio in the polar mesh as shown in [TT] . 




Figure 24: Spherical implosion. Initial polar (left) and unstructured (right) grids. 



Weakly perturbed case with aQ = 2 x 10^^. Now, we investigate the capa- 
bility of our CCALE-MOF algorithm to treat perturbed interfaces on both 
non-structured and polar meshes. To this end, comparisons with pure La- 
grangian results are first achieved for weakly perturbed interfaces imposing 
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Figure 25: Spherical implosion without deformation. Mesh and density for polar (left) 
and unstructured (right) grids at final time tf ~ 3. 



ao = 2 X 10 ^. Here for both polar and hybrid meshes, the GCNS is used. As 



demonstrated on Fig. 26, for the polar mesh as well as for the non-structured 
mesh, ALE results, especially concerning the interface deformation, are in 
very good agreement to thoses obtained thanks to pure Lagrangian compu- 
tations. Furthemore, one should note that for the ALE computation on polar 
grid the quality of the mesh is improved near the origin. Indeed, the central 
cells are not systematically shifted to the origin contrary to computations 
achieved using CNS rezoning. 

Strongly perturbed case with ao = 1 x 10"^. Finally, we perform a computa- 
tion of this implosion for a more pertubated interface choosing five times 
greater than previously with = Ix 10~^. Due to mesh tangling, this is not 
possible to purchase such a test case using only Lagrangian method whose 
computation fails for t > tjau = 2.6. Here, only results obtained thanks to 
our axisymmetric multi-material CCALE-MOF are presented. Contrary, to 



47 



Figure 26: Spherical implosion with small deformation. Mesh and density for Lagrangian 
(top) and ALE (bottom) computations at final time = 3 for both polar (left) and 
unstructured (right) grids. 



Lagrangian computations, the multi-material ALE simulations run without 
any difficulties thanks to specific rezoning. For both grids, final results (see 



Fig. 27) are very close. In particular we note the Ray leigh- Taylor instability 
has grown in a same way leading to similar interface shape deformation at 
final time. 
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Figure 27: Spherical implosion with important deformation. Mesh and density for ALE 
computation for both polar (left) and unstructured (right) grids at final time tf =3. 

7. Conclusion and future work 

In this paper, we have presented several extensions concerning a Cell- 
Centered Arbitrary Lagrangian-Eulerian (CCALE) strategy using the Mo- 
ment of Fluid (MOF) interface reconstruction devoted to the numerical sim- 
ulation of multi-material compressible flows especially in axisymmetric ge- 
ometry on both polar and Cartesian unstructured meshes. To this end, we 
have introduced a simple and unified formulation of the Lagrangian scheme 
relying on an area-weighted formulation, a multi-material MOF interface re- 
construction, a new formulation of rezoning for both polar and Cartesian 
grids and finally a general hybrid remap procedure for both axisymmetric 
and Cartesian geometry. As demonstrated on several academical as well as 
ICF-like test cases, the proposed method remains accurate and robust. 
As future work, we plan to incorporate the proposed method in the multi- 
physic code CHIC dedicated to the simulation of ICF experiment. The main 
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goal is to treat eventually more general configurations notably coupling 
alistic EOS, laser energy deposition, with multi-material hydrodynamics 
the lines of [5]. 
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